In [1]:
%matplotlib inline
In [2]:
# load auxiliary libraries
import numpy as np
import matplotlib.pyplot as plt
from astropy.io import fits
# import stingray
import stingray
plt.style.use('seaborn-talk')
Open the event file with astropy.io.fits
In [3]:
f = fits.open('emr_cleaned.fits')
The time resolution is stored in the header of the first extension under the Keyword TIMEDEL
In [4]:
dt = f[1].header['TIMEDEL']
The collumn TIME
of the first extension stores the time of each event
In [5]:
toa = f[1].data['Time']
Let's create a Lightcurve from the Events time of arrival witha a given time resolution
In [6]:
lc = stingray.Lightcurve.make_lightcurve(toa=toa, dt=dt)
In [7]:
lc.plot()
Let's create a dynamic powerspectrum with the a segment size of 16s and the powers with a "leahy" normalization
In [8]:
dynspec = stingray.DynamicalPowerspectrum(lc=lc, segment_size=16, norm='leahy')
The dyn_ps attribute stores the power matrix, each column corresponds to the powerspectrum of each segment of the light curve
In [9]:
dynspec.dyn_ps
Out[9]:
To plot the DynamicalPowerspectrum matrix, we use the attributes time
and freq
to set the extend of the image axis. have a look at the documentation of matplotlib's imshow()
.
In [10]:
extent = min(dynspec.time), max(dynspec.time), max(dynspec.freq), min(dynspec.freq)
plt.imshow(dynspec.dyn_ps, origin="lower left", aspect="auto", vmin=1.98, vmax=3.0,
interpolation="none", extent=extent)
plt.colorbar()
plt.ylim(700,850)
Out[10]:
In [11]:
print("The dynamical powerspectrun has {} frequency bins and {} time bins".format(len(dynspec.freq), len(dynspec.time)))
In [12]:
print("The current frequency resolution is {}".format(dynspec.df))
Let's rebin to a frequency resolution of 2 Hz and using the average of the power
In [13]:
dynspec.rebin_frequency(df_new=2.0, method="average")
In [14]:
print("The new frequency resolution is {}".format(dynspec.df))
Let's see how the Dynamical Powerspectrum looks now
In [15]:
extent = min(dynspec.time), max(dynspec.time), min(dynspec.freq), max(dynspec.freq)
plt.imshow(dynspec.dyn_ps, origin="lower", aspect="auto", vmin=1.98, vmax=3.0,
interpolation="none", extent=extent)
plt.colorbar()
plt.ylim(500, 1000)
Out[15]:
In [16]:
extent = min(dynspec.time), max(dynspec.time), min(dynspec.freq), max(dynspec.freq)
plt.imshow(dynspec.dyn_ps, origin="lower", aspect="auto", vmin=2.0, vmax=3.0,
interpolation="none", extent=extent)
plt.colorbar()
plt.ylim(700,850)
Out[16]:
Let's try to improve the visualization by rebinnin our matrix in the time axis
In [17]:
print("The current time resolution is {}".format(dynspec.dt))
Let's rebin to a time resolution of 64 s
In [18]:
dynspec.rebin_time(dt_new=64.0, method="average")
In [19]:
print("The new time resolution is {}".format(dynspec.dt))
In [20]:
extent = min(dynspec.time), max(dynspec.time), min(dynspec.freq), max(dynspec.freq)
plt.imshow(dynspec.dyn_ps, origin="lower", aspect="auto", vmin=2.0, vmax=3.0,
interpolation="none", extent=extent)
plt.colorbar()
plt.ylim(700,850)
Out[20]:
Let's use the method trace_maximum()
to find the index of the maximum on each powerspectrum in a certain frequency range. For example, between 755 and 782Hz)
In [21]:
tracing = dynspec.trace_maximum(min_freq=755, max_freq=782)
This is how the trace function looks like
In [22]:
plt.plot(dynspec.time, dynspec.freq[tracing], color='red', alpha=1)
plt.show()
Let's plot it on top of the dynamic spectrum
In [23]:
extent = min(dynspec.time), max(dynspec.time), min(dynspec.freq), max(dynspec.freq)
plt.imshow(dynspec.dyn_ps, origin="lower", aspect="auto", vmin=2.0, vmax=3.0,
interpolation="none", extent=extent, alpha=0.7)
plt.colorbar()
plt.ylim(740,800)
plt.plot(dynspec.time, dynspec.freq[tracing], color='red', lw=3, alpha=1)
plt.show()
The spike at 400 Hz is probably a statistical fluctutations, tracing by the maximum power can be dangerous!
We will implement better methods in the future, stay tunned ;)